# -*- coding: utf-8 -*-
"""
Created on Wed Aug 15 17:12:07 2012

@author: eduardo
"""


    % D-Step: (nD D-Steps for each sample)
    %======================================================================
    for iY = 1:nY
 
 
        % D-Step: until convergence for static systems
        %==================================================================
 
        % derivatives of responses and inputs
        %------------------------------------------------------------------
        pu.z = spm_DEM_embed(Z,n,iY);
        pu.w = spm_DEM_embed(W,n,iY);
        qu.u = spm_DEM_embed(U,n,iY);
        
        % pass action to pu.a
        %------------------------------------------------------------------
        pu.a = qu.a;
 
        % evaluate generative model
        %------------------------------------------------------------------
        [pu dg df] = spm_ADEM_diff(G,pu);
 
        % tensor products for Jacobian
        %------------------------------------------------------------------
        Dgda = kron(spm_speye(n,1,1),dg.da);
        Dgdv = kron(spm_speye(n,n,1),dg.dv);
        Dgdx = kron(spm_speye(n,n,1),dg.dx);
        dfda = kron(spm_speye(n,1,0),df.da);
        dfdv = kron(spm_speye(n,n,0),df.dv);
        dfdx = kron(spm_speye(n,n,0),df.dx);
        
        dgda = kron(spm_speye(n,1,0),dg.da);
        dgdx = kron(spm_speye(n,n,0),dg.dx);
        
        % and pass response to qu.y
        %------------------------------------------------------------------
        for i = 1:n
            y       = spm_unvec(pu.v{i},{G.v});
            qu.y{i} = y{1};
        end
 
        % evaluate recognition model
        %------------------------------------------------------------------       
        [E dE] = spm_DEM_eval(M,qu,qp);
 
        
        % conditional covariance [of states {u}]
        %------------------------------------------------------------------
        qu.c   = spm_inv(dE.du'*iS*dE.du + Pu);
        Hqu.c  = Hqu.c + spm_logdet(qu.c);
        
        % save at qu(t)
        %------------------------------------------------------------------
        qE{iY} = E;
        qC{iY} = qu.c;
        qU(iY) = qu;
        pU(iY) = pu;
 
        % and conditional covariance [of parameters {P}]
        %------------------------------------------------------------------
        ECEu   = dE.du*qu.c*dE.du';
        ECEp   = dE.dp*qp.c*dE.dp';
 
        
        % uncertainty about parameters dWdv, ... ; W = ln(|qp.c|)
        %==================================================================
        if np
            for i = 1:nu
                CJp(:,i)   = spm_vec(qp.c*dE.dpu{i}'*iS);
                dEdpu(:,i) = spm_vec(dE.dpu{i}');
            end
            dWdu  = CJp'*spm_vec(dE.dp');
            dWduu = CJp'*dEdpu;
        end
        
        % change in error w.r.t. action
        %------------------------------------------------------------------
        dyda  = dydv*dgda;
        dydx  = dydv*dgdx;
        for i = 1:n
            Dfdx = kron(spm_speye(n,n,-i),df.dx^(i - 1));
            dyda = dyda + dydx*Dfdx*dfda;
        end
        
        % dE/da with restriction
        %------------------------------------------------------------------
        dE.da = dE.dy*Ra*dyda;
        dE.dv = dE.dy*dydv;
        
        % first-order derivatives
        %------------------------------------------------------------------
        dVdu  = -dE.du'*iS*E - Pu*spm_vec({qu.x{1:n} qu.v{1:d}}) - dWdu/2;
        dVda  = -dE.da'*iS*E - Pa*spm_vec( qu.a{1:1});
        
        % and second-order derivatives
        %----------------------------------------------------------------
        dVduu = -dE.du'*iS*dE.du - Pu - dWduu/2 ;
        dVdaa = -dE.da'*iS*dE.da - Pa;
        dVduv = -dE.du'*iS*dE.dv;
        dVduc = -dE.du'*iS*dE.dc;
        dVdua = -dE.du'*iS*dE.da;
        dVdav = -dE.da'*iS*dE.dv;
        dVdau = -dE.da'*iS*dE.du;
        dVdac = -dE.da'*iS*dE.dc;

        
 
         
        % D-step update: of causes v{i}, and hidden states x(i)
        %==================================================================
 
        % states and conditional modes
        %------------------------------------------------------------------
        p     = {pu.v{1:n} pu.x{1:n} pu.z{1:n} pu.w{1:n}};
        q     = {qu.x{1:n} qu.v{1:d} qu.u{1:d} qu.a{1:1}};
        u     = [p q];  
        
        % gradient
        %------------------------------------------------------------------
        dFdu  = [                              Dp*spm_vec(p); 
                 spm_vec({dVdu; dVdc; dVda}) + Dq*spm_vec(q)];
 
 
        % Jacobian (variational flow)
        %------------------------------------------------------------------
        dFduu = spm_cat(...
                {Dgdv  Dgdx Dv   []   []       []    Dgda;
                 dfdv  dfdx []   dfdw []       []    dfda;
                 []    []   Dv   []   []       []    [];
                 []    []   []   Dx   []       []    [];
                 dVduv []   []   []   Du+dVduu dVduc dVdua;
                 []    []   []   []   []       Dc    []
                 dVdav []   []   []   dVdau    dVdac dVdaa});
 
 
        % update states q = {x,v,z,w} and conditional modes
        %==================================================================
        du    = spm_dx(dFduu,dFdu,td);
        u     = spm_unvec(spm_vec(u) + du,u);
 
        % and save them
        %------------------------------------------------------------------
        pu.v(1:n) = u([1:n]);
        pu.x(1:n) = u([1:n] + n);
        qu.x(1:n) = u([1:n] + n + n + n + n);
        qu.v(1:d) = u([1:d] + n + n + n + n + n);
        qu.a(1:1) = u([1:1] + n + n + n + n + n + d + d);

    end % sequence (nY)